function res = trainModel_libsvm4_slle(stego, cover, ntrain, ntest, loop, settings)
% TRAINMODEL train SVM model using LIBSVM.
% Usage:
%   RES = TRAINMODEL_LIBSVM4_SLLE(STEGO, COVER, NTRAIN, NTEST, LOOP, SETTINGS)
%
%  This function splits training sets and testing sets, trains the model, 
%  and return the testing result. 
%  Selected stego and cover images are in pairs.
%
%  The following is equivalence:
%   RES = TRAINMODEL_LIBSVM4_SLLE(STEGO, COVER, NTRAIN, NTEST)
%   RES = TRAINMODEL_LIBSVM4_SLLE(STEGO, COVER, NTRAIN, NTEST, 10)
%   RES = TRAINMODEL_LIBSVM4_SLLE(STEGO, COVER, NTRAIN, NTEST, 10, SETTINGS)
%  where SETTINGS is
%   s = struct('verbose', 0, 'c', nan, 'g', nan, 'ncg', 0, ...
%   'kerneltype', 0, 'svmtrain','svmtrain','svmpredict','svmpredict')
%
% Input:
%   STEGO  : a L*N matrix for features of stego images.
%   COVER  : a M*N matrix for features of cover images.
%   NTRAIN : size of training set. (one pair of stego/cover images is 2!)
%   NTEST  : size of testing set.
%   LOOP   : number of loops. Default value is 10.
%   SETTINGS: a struct of various settings :
%     verbose    : 3 show debug messages.
%                  2 show results, errors and steps.
%                  1 show errors and steps.
%                  0 only show error messages.
%     c, g       : Value for LIBSVM. Ignore it, value 'automatic' or nan
%            tells the program to find  best values for C,G using search
%            (and this process is very slow).
%     ncg        : size of training sample when searching for C, G. If the
%            value is 0, it is defined as NTRAIN.
%     kerneltype : See kernel type for svmtrain. Default value is 0
%             (linear).
%     svmtrain   : path to svmtrain executable.
%     svmpredict : path to svmpredict exectuable.
%     title      : some name.
%     parallel   : run in parallel (for example, in CG selecting, training)
%     KK, DD, a  : for slle.
%  Reference to SVMCG3 for details.
%     vfold      : v-fold for training C and G.
%     cbin       : range of C values.
%     gbin       : range of Gamma values.
% 
% Output:
%   RES : A struct with the following fields:
%       CG_OPT      : Option used in C, G selecting.
%       TRAIN_OPT   : Option used in trainging.
%       PREDICT_OPT : Option used in predicting.
%       
%       MEAN_ERR    : Mean probability of predicting errors.
%       MEAN_FA     : Mean of false alarm. (0 predicted as 1)
%       MEAN_MD     : Mean of 
%       MIN_ERR     : Minimal error among all loops.
%
% See also:
%   SVMCG3
%   TRAINMODEL_LIBSVM2

%% parse inputs
if nargin < 6
    settings = struct();
end

if ~isfield(settings, 'title')
    settings.title = '_';
end
if ~isfield(settings, 'svmtrain')
    settings.svmtrain = 'svm-train';
end
if ~isfield(settings, 'svmpredict')
    settings.svmpredict = 'svm-predict';
end
if ~isfield(settings, 'verbose')
    settings.verbose = 0;
end
if ~isfield(settings, 'kerneltype')
    settings.kerneltype = 0;
end
if ~isfield(settings, 'ncg')
    settings.ncg = 0;
end
if ~isfield(settings, 'c') || ~isfield(settings, 'g')
    settings.c = nan;
    settings.g = nan;
end
if ~isfield(settings, 'pair')
    settings.pair = 1;
end

if ~isfield(settings, 'parallel')
    settings.parallel = 0; % disable it by default, experimental
end

if ~isfield(settings, 'a')
    settings.a = 1;
end
settings.cg_opt = ' -h 1 -m 512';
settings.train_opt = '-h 1 -m 512';
settings.predict_opt = '';

settings.cg_opt = ['-t ',num2str(settings.kerneltype),' ', settings.cg_opt];
settings.train_opt = ['-t ',num2str(settings.kerneltype),' ', settings.train_opt];
if settings.verbose == 0
    settings.train_opt = ['-q ', settings.train_opt];
%    settings.predict_opt = ['-q ', settings.predict_opt];
    settings.predict_opt = [settings.predict_opt];
end

if settings.ncg == 0
    settings.ncg = ntrain;
end

%% load feature matrix
if ischar(stego)
    if settings.verbose >= 2
        fprintf('loading stego feature from file %s.\n', stego);
    end
    load(stego,'-mat','res');
    stego = res;
    clear res;
end
[L,N] = size(stego);
if ischar(cover)
    if settings.verbose >= 2
        fprintf('loading cover feature from file %s.\n', cover);
    end
    load(cover,'-mat','res');
    cover = res;
    clear res;
end
[M,N1] = size(cover);
if (N ~= N1)
    error('trainModel:featureNumNotMatch','Number of features in STEGO and COVER does not match.');
end
clear N1
if (L ~= M) && (settings.pair)
    error('Number of samples in stego and cover are not equal, set SETTINGS.pair to 0. This may cause suboptimal classification results.');
    settings.pair = 0;
end

if nargin < 5
    loop = 10;
elseif loop <= 0
    error('trainModel:invalidInput','value of loop must be >0.');
end

%% initialize output value
res = struct('CG_OPT', settings.cg_opt, 'TRAIN_OPT', settings.train_opt, ...
    'PREDICT_OPT', settings.predict_opt, 'STEGO_SAMPLE_NUM', M, ...
    'COVER_SAMPLE_NUM', L, 'FEATURE_DIM', N, 'USE_PAIR', settings.pair, ...
    'PARALLEL', settings.parallel);

nSample = L+M;
% stego(1) | cover(0)
if strcmpi(settings.c, 'automatic') || strcmpi(settings.g, 'automatic') ...
        || isnan(settings.c) || isnan(settings.g)
    fprintf('finding best c and g\n');
    nt = floor(settings.ncg / 2);
    % nt from L and M
    p1 = randperm(L);
    if settings.pair
        p2 = p1;
    else
        p2 = randperm(M);
    end
    train = [stego(p1(1:nt), :); cover(p2(1:nt), :)];
    ltrain = [ones(nt,1); zeros(nt,1)];
    train = slle(train',ltrain,settings.KK,settings.DD,settings.a);
    if settings.verbose >= 2
        fprintf('train set: stego = cover = %d\n', nt);
    end
    tic;
    [acc, c, g, cvmat] = svmcg3(ltrain, train, settings);
    if settings.verbose
        disp(cvmat);
        fprintf('c=%g, g=%g\n', c, g);
    end
    T=toc;
    settings.c = c;
    settings.g = g;
    [res.CG_ACC, res.CG_TIME, res.CG_MAT] = deal(acc, T, cvmat);
    clear train ltrain
end

if (ntrain <= 0) || (ntest < 0) || (ntrain + ntest > nSample)
    error('trainModel:invalidInput', 'NTRAIN + NTEST must not be bigger than total sample numbers.');
end
[res.NUM_TRAIN, res.NUM_TEST, res.NUM_LOOP] = deal(ntrain, ntest, loop);
[res.C, res.G] = deal(settings.c, settings.g);
[res.MEAN_ERR, res.MEAN_FA, res.MEAN_MD, res.MIN_ERR] = deal(nan);

errors_fa = nan(loop, 1); %false alarm 0 said to be 1
errors_md = nan(loop, 1); %miss detect 1 said to be 0
accs = nan(loop, 3);
errors = nan(loop, 1);

ntrain = floor(ntrain / 2);
ntest = floor(ntest/2);
output = '';

if settings.verbose
    fprintf('trn/tst = %d/%d, loop = %d\n',ntrain, ntest, loop);
    if settings.pair
        fprintf('train set: %d pairs\n', ntrain);
        fprintf('test set : %d pairs\n', ntest);
    else
        fprintf('train set: %d stego/cover images\n', ntrain);
        fprintf('test set : %d stego/cover images\n', ntest);
    end
end

%% train and test
% implementation of LIBSVM is not thread safe
es = 0;
tic;
for i = 1:loop
    p1 = randperm(L);
    if settings.pair
        p2 = p1;
    else
        p2 = randperm(M);
    end
    
    if settings.verbose
        fprintf('Round %d/%d\n',i,loop);
    end
    % 1 .. ntrain              => train
    % ntrain+1 .. ntrain+ntest => test
    % ntrain+ntest+1 ...       => skipped
    train = [stego(p1(1:ntrain), :); cover(p2(1:ntrain), :)];
    ltrain = [ones(ntrain,1);zeros(ntrain,1)];
    old_train = train;
    train = slle(train',ltrain,settings.KK,settings.DD,settings.a)';

    t1=cputime;
    model = libsvmtrain_f(ltrain, train, settings, settings.c, settings.g, 1);
    t2=cputime-t1;
    fprintf('time=%g\n',t2);
    
    test = [stego(p1(ntrain+1:ntrain+ntest), :); cover(p2(ntrain+1:ntrain+ntest), :)];
    ltest = [ones(ntest,1);zeros(ntest,1)];
    test = deTestD(old_train',test',train', settings.KK)';

    [tr, acc, ~] = libsvmpredict_f(ltest, test, model, settings, 1);
    clear test train ltrain old_train

    pfa = sum((tr == 1) & (ltest == 0)) / sum(ltest == 0);
    pmd = sum((tr == 0) & (ltest == 1)) / sum(ltest == 1);
    err = (pmd + pfa)/2;
    errors_fa(i) = pfa;
    errors_md(i) = pmd;
    errors(i) = err;
    es = es + err;
    accs(i,:) = acc';
    %if settings.verbose < 1
    %    fprintf(repmat('\b',1,length(output)));
    %end
    output = sprintf('%d/%d : acc = %0.4f / %0.4f\n', i, loop, 1-err, 1-es/i);
    fprintf(output);
end
fprintf('\n');
T = toc;
res.TIME = T;
res.errors_fa = errors_fa;
res.errors_md = errors_md;
res.errors = errors;
res.accuracy = accs;
res.MIN_ERR = min(errors);
errors(isnan(errors)) = 0;
errors_fa(isnan(errors_fa)) = 0;
errors_md(isnan(errors_md)) = 0;
[res.MEAN_ERR, res.MEAN_FA, res.MEAN_MD] = deal(sum(errors) / loop, sum(errors_fa) / loop, sum(errors_md) / loop);
if settings.verbose
    disp(res)
    fprintf('Mean Accuracy: %0.2f%%\n',(1-res.MEAN_ERR)*100);
end 
end